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EXECUTIVE  SUMMARY 


OBJECTIVE 

The  objective  is  to  provide  environmental  support  to  system  designers  and  field 
commanders  to  accurately  forecast  pertinent  conditions  and  to  make  real  time  assessment/ 
tactical  conditions. 

RESULTS 

A  raytrace  technique  has  been  developed  to  show  the  wave  front  path  as  it  propagates 
through  the  laterally  heterogeneous  medium  where  the  index  of  refraction  is  allowed  to  vary 
both  vertically  and  horizontally. 

RECOMMENDATION 

This  effort  is  specially  geared  toward  future  Navy  systems  that  are  presently  in  the 
planning  stage  or  under  development. 
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INTRODUCTION 


Radio  wa\es  travel  in  straight  lines  through  any  isotropic,  homogeneous  medium  with 
a  velocity  which  depends  upon  the  speed  of  light  in  a  vacuum  and  the  refractive  index  of  the 
medium.  If  the  index  of  refraction  changes,  the  velocity  and  possibly  the  direction  of  travel 
will  change. 

Over  major  land  and  ocean  surfaces,  conditions  of  temperature  and  moisture  are 
sufficiently  homogeneous  to  cause  the  formation  of  air  masses  which  provide  a  laterally 
homogeneous  medium  for  electromagnetic  wave  propagation.  Since  air  mass  conditions  are 
the  norm  for  most  of  the  troposphere,  many  raytrace  techniques  have  been  developed  to 
illustrate  a  wave  front  as  it  propagates  through  this  laterally  homogeneous  medium. 

At  air  mass  boundaries  associated  with  wave  cyclones.  Foehn  circulations,  land  ocean 
interfaces,  etc  .  the  conditions  of  lateral  homeogeneity  often  do  not  exist.  A  raytrace 
technique  has  been  developed  to  show  the  wave  front  path  as  it  propagates  through  this 
laterally  heterogeneous  medium  where  the  index  of  refraction  is  allowed  to  vary  both 
vertically  and  horizontally.  Appendix  A  lists  the  program  code  (Microsoft  QuickBASIC 
version  2.0 1 )  for  the  raytrace  technique  as  employed  upon  an  IBM  PC.  XT,  AT  or  compatible 
using  the  MS-DOS  version  3.2  operating  system. 


BACKGROUND 

I  he  retractive  index  n  of  a  parcel  of  air  is  defined  as  the  ratio  of  the  velocity  of 
propagation  of  an  electromagnetic  (EM)  wave  in  a  vacuum  to  that  in  the  air.  Since  the 
retractive  index  of  the  atmosphere  is  slightly  greater  than  unity.  EM  waves  travel  slightly 
slower  in  air  than  in  a  vacuum.  Close  to  the  earth's  surface,  the  numeric  value  of  the  refractive 
index  is  usually  between  1.00025  and  1.0004.  For  convenience,  the  refractivity  V  is  defined  by 
Bean  and  Dutton  i  I96M  as 

V  =  (/j  I )  I()6  (1) 

such  that  surface  values  ot  refractivity  vary  between  250  and  400.  Refractivity  may  be 
expressed  as  a  function  ot  atmospheric  pressure,  temperature,  and  humidity  by  the  relationship 

7Tb  />  3,73  l()5e 

- *  -  (2) 

r 


vv  here 

P  -  atmospheric  pressure  (millibars) 

/  -  atmospheric  temperature  (  Kelv  in) 

(’  -  atmospheric  water  vapor  pressure  (millibars)  . 

•\  ray  path  under  well-mixed  (standard)  atmospheric  conditions  will  bend  downward 
ai  a  late  less  than  the  curvature  ot  the  earth.  Certain  atmospheric  conditions  could  also  lead 
In  the  lav  path  bending  downward  at  a  rate  exceeding  the  curvature  ot  the  earth  (trapping), 
downward  moic  than  standard  but  not  sutticient  tor  trapping  (super-retraction),  or  upward 
(sub-rcliactionl  figure  I  illustrates  these  retractive  conditions. 


I 


Figure  1  Relative  bending  for  four  types  of  refractive  conditions. 


As  a  further  convenience  in  computing  and  graphically  displaying  these  types  of 
refractive  conditions,  a  modified  refractivity  M  may  be  defined  as 


M  = 


(" 


1  +  a)I()6 


(3) 


where 

a  -  mean  earth's  radius  and 
r  =  specified  height  above  the  earth's  surface. 

The  advantage  of  the  modified  refractivity  is  that  trapping  environments  are  indicated 
with  a  negative  M-unit  gradient  with  increasing  height.  Table  I  lists  the  relations  between  ,V, 
M,  and  the  refractive  conditions. 


Table  I.  The  relationship  between  ,Y,  M,  and  refractive  conditions. 


Types  of  refraction 


N-unit  gradient  (N  km) 


M-unit  gradient  (M  km) 


Trapping 

Super-refraction 

Standard 

Sub-refraction 


<  1 57 
157  to  79 

79  to  0 

>  0 


<  0 
0  to  79 


RAY  THEORY 


The  propagation  of  EM  waves  within  the  troposphere  may  be  simply  explained  by  ray 
theory  as  developed  by  Reed  and  Russell  ( 1966).  As  state above,  the  refractive  conditions 
lead  to  a  change  in  the  propagating  ray's  direction.  Figure  2  illustrates  the  direction  change 
and  is  defined  by  Snell’s  law  : 

n  |  cos(a  | )  =  ri2  cos(a->)  (4) 

where  the  boundary  is  a  plane  surface. 


Figure  2.  Refraction  at  a  boundary 


In  the  atmosphere,  the  refractive  index  changes  continuously  or  may  be  thought  to 
change  linearly  when  the  atmosphere  is  divided  into  an  infinite  number  of  parallel  boundaries 
and  the  refractive  index  changes  by  infinitesimal  amounts.  As  illustrated  in  figure  3.  Snell's 
law  may  now  be  written  as 

n |  costa | )  =  cos(a(j)  (5) 

where  /ij  and  a|  are  functions  of  height  and  Hq  and  oq  are  fixed  values  at  a  given  reference. 


Figure  3  Refraction  through  an  infinitesimally  changing  n 


Since  the  earth's  surface  is  spherical  and  not  planar,  Snell's  law  may  be  rewritten  as 


rt|  rj  costa | )  =  rig  rg  cos(oq) 


(6) 


where 

r  |  and  rg  are  the  distances  from  the  earth's  center  to  the  boundaries  and  aj  and  ag 
are  the  angles  between  the  ray  and  the  planes  normal  to  the  radius  vector  at  the  points  where 
the  ray  crosses  the  boundaries  as  shown  by  figure  4. 


Figure  4.  Illustration  of  Snell’s  law  for  a 
curved  earth. 

If  z  is  the  height  of  an  antenna  above  the  earth's  surface  and  a  is  the  radius  of  the 
earth,  then 

r |  -  z  *  a  and  rg  =  a  (7) 

and  equation  6  mav  be  written  as 

rt  |  (  I  +  'jj’j  cos(o|)  =  rjg  cos(ftg) .  (8) 

W  ith  the  assumption  that  n  is  close  to  units  and is  a  vary  small  quantity  compared  to  unity 

,9' 

and  lor  small  values  of  « 

■> 

cost  a)  =  I  (10) 
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'mm 


Snell's  law  may  be  written  as 


l(al  Qi)  =  ”!  "0+^  '  <H) 

Substituting  the  modified  refractivity  lor  the  refractive  index,  solving  for  the  ending 
height  and  angle  of  a  ray  when  given  the  beginning  position,  and  the  ray's  initial  launch  angle 
yields  the  traditional  raytrace  equations  employed  by  well-established  rayiracing  theory. 

These  equation  are 


and 


A  V/  , 

°i  =  fto  +  —  (vl 


V  10  6 


(12) 


I 


■0  A_U 
Ar 


10 


(13) 


w  here 


ftp  =  ray  angle  at  start  of  calculation  (radians) 
o  I  =  ray  angle  at  end  of  calculation  (radians) 

Cq  =  height  at  start  of  calculation  (meters) 

C|  -  height  at  end  of  calculation  (meters) 
x'O  -  range  at  start  of  calculation  (kilometers) 

V|  -  range  at  end  of  calculation  (kilometers) 

-  change  in  M-units  with  respect  to  a  change  in  height  (M  units  per  meter). 


LINEAR  MODEL  OF  A  HETEROGENEOUS  ATMOSPHERE 

I  o  model  an  atmosphere  which  varies  both  vertically  and  horizontally,  it  is  necessary 
to  make  assumptions  concerning  the  behav  ior  of  the  refractive  index.  At  air  mass  boundaries 
caused  by  large  scale  subsidence,  such  as  within  the  trade  wind  inversion  as  illustrated  bv 
Hitney.  et  al.  ( 1985)  in  figure  5.  it  is  natural  to  assume  a  horizontally  homogeneous  continua¬ 
tion  of  the  boundary  With  vertical  air  mass  boundaries,  such  as  those  associated  with  mid- 
latitude  wave  cyclones  as  shown  by  Bean  and  Dutton  ( 1968)  in  figure  6.  or  those  associated 
with  radiational  cooling  as  shown  by  Morrissey  (1985)  in  figure  7.  the  horizontal  variation  of 
reft  activity  may  be  more  complicated 

In  this  raytrace  technique,  it  is  assumed  the  atmosphere  can  be  divided  into  layers 
with  piece-wise  linear  variations  in  the  layer  height  and  piece-wise  linear  refractivity  within 
each  layer  In  areas  where  a  duct  degenerates  to  a  linear  profile  or  a  duct  splits  to  form  several 
detached  ducts,  layer  boundaries  a  re  determined  by  the  meteorological  conditions. 
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20  40  60 


Figure  5.  Space  cross  section  in  B-units,  12  March  1948,  San  Diego,  CA  seaward  to 


Guadalupe  Islands.  ^  B  -  £(«  -  I )  +  I06.^ 
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Figure  6.  Space  cross  section  in 
N  units,  1  500Z,  1  9  February  1  952 
Solid  vertical  line  represents 
position  of  mid-latitude  wave 
cyclone  cold  front. 
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Knowing  the  vertical  structure  of  refractivity  and  the  horizontal  distance  between  two 

A  \1 

refractivity  profiles  "a"  and  “b,”  the  layer  boundary  height  r.  and  the  refractivity  gradient 
is  given  by  the  relationship 


(~6  ) 

(xb  -  '  a  ) 


(ft»  ft») 


A  A/  ,  , 
v+— 
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OPERATOR  /  PROGRAM  INTERACTIONS 


ut\ 
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DATA  REQUIREMENTS 

The  operator  is  allowed  to  interact  with  the  program  in  two  ways.  First,  the  operator 
may  enter  data  from  the  keyboard  or  request  the  data  be  read  from  a  disk  file  previously 
created.  The  input  data  are 

a.  environmental  data  consisting  of  triplets  of  pressure,  temperature,  and  relative 
humidity,  pairs  of  N-units  and  height,  or  pairs  of  M-units  and  height. 

b.  transmitter  height 

c.  radiating  antenna  beamwidth 

d.  .  radiating  antenna  elevation  angle 

e.  the  number  of  rays  to  draw 

Second,  the  operator  may  examine  the  refractivity  profiles  to  manually  determine  the 
layer  boundaries. 


ENVIRONMENTAL  DATA  CONSIDERATIONS 

For  this  program,  the  maximum  number  of  environmental  profiles  which  may  be 
entered  is  five.  This  number  was  chosen  for  graphics  convenience  only.  This  raytracc 
technique  may  be  applied  to  any  number  of  profiles  needed  to  fully  describe  the  laterally 
heterogeneous  atmosphere. 

Prior  to  the  determination  of  layer  boundaries,  the  environmental  data  profiles  must 
be  modified  to  insure  each  profile  starts  at  a  height  of  zero  and  each  profile  terminates  at  the 
same  elevation.  This  procedure  is  accomplished  by  the  following  steps: 


a.  If  the  environmental  data  is  entered  as  triplets  of  pressure,  temperature,  and 
relative  humidity,  convert  the  data  to  M-units  (equations  2  and  3),  and  height  using  the 
relationship  derived  by  Berry  (1945): 

p  =  n0  (  i  (16) 

where 

P  =  atmospheric  pressure  at  a  height  of  z 
P  -  atmospheric  pressure  at  a  height  of  zero 
/3  =  atmospheric  lapse  rate 

T  -  temperature  at  a  height  of  zero 
g  =  gravity 

R  =  moist  atmospheric  gas  constant 
c  =  height, 

b.  If  the  environmental  data  is  entered  as  pairs  of  N-units  and  height  convert  the 
data  to  M-units  (equation  3). 

c.  If  necessary,  the  lowest  height  within  each  profile  is  set  to  zero  and  a  surface 
M-unit  value  is  extrapolated  using  the  standard  atmospheric  M-unit  gradient  of  1 18  M-units 
per  kilometer  as  specified  by  Bean  and  Dutton  ( 1968). 

d.  If  necessary,  extend  the  profile  height  to  the  height  of  the  highest  profile  (/max) 
and  determine  a  corresponding  M-unit  \alue  from  the  exponential  atmosphere  model  of  Bean 
and  Dutton  ( 1968): 


M, 

“max 


=  157  +  \f 


sfc  e 


(17) 


w  here 

Ms)c  =  M-unit  value  at  a  height  of  zero 
M^m  =  M-unit  \alue  at  a  height  of  I  kilometer. 

e.  An  M-unit  gradient  is  determined  at  each  profile  height  using  the  relationship 


a.v/  (A/.+ 1  v/') 

A:i  (Ji+I  ri) 


(18) 


8 


HEIGHT 


where  the  M-unit  gradient  at  the  profile  top  (7max)  is  given  by  differentiating  equation  17 
with  respect  to  height 


LAYER  BOUNDARY  DETERMINATION 

The  modified  refractivity  versus  height  profiles  are  computed  and  displaced  on  the 
terminal  screen  as  illustrated  in  figure  8.  Through  screen  prompts  and  the  use  of  a  moveable 
cursor,  the  operator  is  allowed  to  specify  the  layer  boundaries.  Each  profile  must  contain  an 
equal  number  of  layers.  Figure  9  illustrates  the  operator  specified  boundaries  for  the 
complicated  profile  set  of  figure  8.  Since  the  layer  boundaries  pl^v  a  major  role  in  this 
raytrace  method,  the  operator  must  insure  the  boundaries  represent  the  most  reasonable  and 
meteorologically  correct  assessment.  As  demonstrated  in  figure  8.  the  elevated  and 
surface-based  trapping  layers  extend  throughout  all  the  profiles  allowing  for  a  simple 
assessment.  The  mid-level  trapping  layer  however,  is  absent  from  profiles  I  and  2.  Since  the 
trapping  layer  was  decreasing  in  thickness  from  profiles  5  to  3.  the  operator  chose  to  eliminate 
the  layer  by  allowing  the  M-unit  gradient  to  change  from  standard  to  trapping  between 
profiles  2  and  3.  Several  methods  of  automating  this  boundary  determination  process, 
including  the  use  of  potential  temperature,  are  the  subject  of  further  investigation. 


M  UNIT 
PROFILE  1 


M  UNIT 


M  UNIT 


M  UNIT 
PROFILE  4 


M  UNIT 
PROFILE  5 


Figure  8  An  example  heterogeneous  atmosphere  defined  by  5  M-unit  versus  height  profiles 
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M-UNIT  M  UNIT  M  UNIT  M  UNIT  M-UNIT 


Figure  9  Example  layer  boundaries 


RAYTRACE  TECHNIQUE 

The  following  steps  are  taken  in  the  execution  of  the  raytrace: 

a.  Establish  a  beginning  range  (vg)  of  zero  and  a  beginning  height  (eg)  equal  to  the 
transmitter  height.  The  first  ray's  beginning  launch  angle  (ag)  is  set  to  half  of  the  beamwidth 
below  the  antenna  elevation  angle. 

b.  Compute  a  M-unit  gradient  at  this  beginning  point  by  using  equation  15. 

c.  I  sing  equations  12  and  1 3,  compute  an  ending  eight  (C|)  and  angle  (aj ).  by 
using  as  an  ending  range  (,V| ).  the  starting  range  plus  an  increment  equal  to  I  25th  of  the 
total  range  between  the  profiles  being  considered.  The  increment  of  range  was  chosen  to 
produce  a  smooth  curve  appearance  of  the  ray  upon  the  graphical  program  output.  In 
subsequent  range  steps,  the  ending  range  must  be  compared  to  the  total  range  between 
profiles.  If  the  ending  range  exceeds  the  range  between  profiles,  the  ending  range  is  set  to  this 
total  range  with  a  new  ending  height  and  angle  being  computed. 

d.  The  ending  angle  is  examined.  If  it  shows  a  sign  reversal  from  the  starting  angle, 
the  ray  has  passed  through  a  maximum  or  minimum  point.  Should  this  be  the  case,  the 
ending  angle  is  set  to  zero  and  a  range  and  height  of  this  maximum  or  minimum  is  computed 
using  equations  12  and  13.  As  started  in  (c)  above,  the  ending  range  must  be  compared  to  the 
total  range  between  the  profiles  being  considered  and  any  necessary  adjustments  made. 

e.  At  the  ending  range,  the  layer's  upper  and  lower  boundary  heights  are  computed 
and  compared  to  the  ending  height.  If  the  ending  height  is  outside  the  layer,  the  range  (.Vj) 
and  the  height  (rj)  of  the  ray  boundary  interception  are  computed  with  the  relationship 


where 


A.v  = 


a 

$ 

$ 


A  r  =  (</K.Vq  +  A.v)  +  y(  c0 


(21) 

(22) 

(23) 

(24) 


v  -  /0  <t>xQ  y  (25) 

<b  =  slope  of  the  boundary 
y  -  height  of  the  boundary  at  the  profile  range. 

Figure  10  illustrates  the  geometry  and  definition  of  terms  at  the  interception  point. 
Since  the  solution  for  A.v  consists  of  two  possibilities,  the  appropriate  choice  is  the  A.v  which 
is  positive  and  less  than  V|. 


PROFILE  a  PROFILE  b 

Figure  10  Geometry  for  a  ray  crossing  an  environmental  layer's  boundary 

f  I  he  ending  height  is  examined  If  it  is  zero,  the  ray  has  reached  the  ground.  In 
this  case,  the  sign  of  the  ending  angle  is  reversed  to  indicate  a  surface  reflection  of  the  ray . 

g  The  beginning  angle,  height,  and  range  are  reinitiali/ied  with  the  ending  angle, 
height,  and  range  and  steps  h  through  g  are  repeated.  I  his  process  continues  until  either  the 


ray  has  reached  the  maximum  range  between  the  first  and  last  profile  or  has  reached  the 
maximum  height  for  the  environmental  profiles.  At  that  time,  step  a  is  repeated  for  the  next 
ray.  increasing  the  initial  ray's  launch  angle  with  an  increment  equal  to  the  total  beamwidth 
divided  by  the  number  of  rays  to  be  drawn. 


RAYTRACE  EXAMPLES 

Figure  1 1  illustrates  the  raytrace  technique  under  a  homogeneous  environment  with  a 
single  elevated  trapping  layer  between  1000  and  1800  feet.  Table  2  lists  the  transmitter  and 
environmental  data.  This  example  is  equivalent  to  one  derived  from  any  standard  raytrace 
technique. 


RAYTRACE  —  HETEROGENEOUS  ATMOSPHERE 


RANGE  (nmi) 

Figure  1 1  Ray  family  for  laterally  homogeneous  trapping  layer 


Table  2.  Environmental  and 


transmitter  data  used  in  figure  1 1. 


Figures  12  and  13  illustrate  the  ray  trace  tor  a  transmitter  located  w  it  hin  a  single 
trapping  laver.  where  the  layer's  elevation  is  allowed  to  rise  and  tall,  respectively,  with  range. 
The  transmitter  and  environmental  data  are  listed  in  tables  3  and  4  For  comparision 
purposes,  these  data  are  identical  to  that  used  bv  Gurnard,  et  al  ( I9b5).  in  their  work  with 
rising  and  tailing  ducts  and  as  illustrated  in  figures  14  and  15.  I  here  exists  good  agreement 
between  the  traces  obtained  from  this  technique  and  that  described  bv  (iuinard  I  he  ma|or 
advantage  of  this  raytracc  technique  over  that  described  bv  (iuinard  is  there  is  no  need  for  a 
secondary  rav  trace  technique  near  a  "caustic  ”  A  caustic  is  defined  as  the  point  at  which  the 
rays  undergo  a  focusing  and  the  assumptions  for  rav  optics  theory  arc  not  valid 

RAYTRACE  HETEROGENEOUS  ATMOSPHERE 
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Figure  1  2  Ray  family  for  a  rising  trapping  layer 
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Figure  13  Ray  family  Tor  a  falling  flapping  layer 


***;•* 


Table  3.  Environmental  and  transmitter  data  used 
in  figure  12. 


y,y 


Profile  1 

Profile  2 

Height  (ft)  M-units 

Height  (ft) 

M-uniS 

0.0 

350.0 

0.0 

350.0 

3500.0 

483.0 

5500.0 

559.0 

4000.0 

457.0 

6000.0 

553.0 

10000.0 

685.0 

10000.0 

685.0 

Transmitter  height  =  3800  ft 

Antenna  beamwidth  =  1.0° 

Antenna  elevation  angle  =  0° 

Table  4.  Environmental  and  transmitter  data  used 
in  figure  1 3. 


Profile  1 

Profile  2 

Height  (ft)  M-untts 

Height  (ft) 

M-units 

0.0 

350.0 

0.0 

350.0 

3500.0 

483.0 

1500.0 

407.0 

4000.0 

457  0 

2000.0 

381.0 

10000  0 

685.0 

10000  0 

685.0 

I  ransmitter  height 

=  38(H)  It 

Antenna  bcamwidth  =  I  t) 

Antenna  elev ation 

angle  -  0 

V 

.\v; 


HAN<jF  n'"i 


f  iqur.-  1  b  f.imiK  < , f  r,iv  paths  <  ali.ulated  for  a  bOO  ft  thick  interface 
with  ,i  slop*'  ul  1  83  «  10  3  ,t  gradient  m  th,-  interface  of  0  1 
N  un.ts  pcf  ft  anil  ,i  qr  adiont  above  and  below  tin*  interface  of  0  01 

N  unit-,  tier  f! 


Figures  16  and  17  illustrate  the  raytrace  using  actual  measured  atmospheric  conditions 
as  shown  in  figure  5.  For  simplicity  of  demonstration,  the  measured  M-unit  versus  height  data 
were  reduced  to  a  set  of  simple  trilinear  profiles.  Table  5  lists  these  data.  It  can  be  seen  that 
when  the  transmitter  is  located  at  a  boundary  (500  feet  in  figure  16).  the  raytrace  technique  is 
able  to  function  w  ithout  special  considerations  when  the  ray  launch  angle  and  boundary  slope 
are  parallel  to  the  earth's  surface  and  equal  to  zero. 
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Figure  1  6  Ray  family  for  environment  of  figure  5  Transmitter  height  of  500  ft 
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Figure  1  7  Ray  family  for  environment  of  figure  5  Transmitter  height  of  1 000  ft 


Table  5.  Environmental  and  transmitter  data  used 
in  figures  16  and  1 7. 


Profile  1 

Profile  2 

Height  (ft) 

M-units 

Height  (ft) 

Vf-units 

0.0 

335.0 

0.0 

335.0 

500  0 

353.0 

500  0 

353.0 

1 100.0 

327.0 

1 100.0 

327.0 

3000.0 

4  ISO 

3000.0 

418.0 

Profile  3 

Profile  4 

Height  (ft) 

M-units 

Height  (ft) 

M-units 

0.0 

335.0 

0.0 

335.0 

700.0 

360.0 

1  150.0 

376.0 

1500.0 

346  0 

2000.0 

364  0 

3000.0 

4 18.0 

3000.0 

412.0 

Profile  5 

Height  (ft) 

M-units 

0.0 

335.0 

1 150.0 

376.0 

2150.0 

370.0 

3000.0 

4110 

Transmitter  height  =  500  and  1000  ft 

Antenna  beamwidth  -  0.5' 

Antenna  elevation  angle  =  0C 

Distance  between  profiles  =  40  nmi 

A  limitation  to  the  practical  application  of  any  heterogeneous  atmosphere  raytracing 
technique  lays  in  obtaining  sufficient  environmental  data  to  adequately  describe  the  refractive 
conditions.  For  this  reason,  a  simulated  environment,  as  shown  by  figure  8  and  as  listed  in 
table  6.  is  used  to  demonstrate  the  ability  of  the  raytrace  technique  to  handle  a  complicated 
environment.  Figure  18  represents  a  transmitter  located  within  a  surface-based  trapping  layer 
which  decrease  in  thickness  with  range.  It  may  be  seen  that  an  increasing  penetration  angle 
compared  to  a  decreasing  boundary  height  gives  rise  to  energy  "leaking"  from  the  trapping 
layer.  In  figure  19.  the  distance  between  the  fourth  and  fifth  profiles  has  been  increased  to 
more  vividly  demonstrate  that  a  ray  lanched  within  an  elevated  trapping  layer  is  capable  of 
leaving  the  layer  and  being  trapped  within  another  layer. 


Table  6.  Environmental  and  transmitter  data  used 
in  figures  18  and  19. 


Profile  1 

Profile  2 

Height  (ft) 

M-units 

Height  (ft) 

M-units 

0.0 

335.0 

0.0 

335.0 

1400.0 

321.0 

1200.0 

323.0 

3500.0 

396.0 

3500.0 

405.0 

3800.0 

407  0 

3800.0 

416.0 

6000.0 

486.0 

7000.0 

531.0 

7300.0 

473.0 

7800.0 

523.0 

8000.0 

498.0 

8000.0 

530.0 

Profile  3 

Profile  4 

Height  (ft) 

M-units 

Height  (ft) 

M-units 

0.0 

335  0 

0.0 

335.0 

1000  0 

325.0 

700.0 

328.0 

3500.0 

415.0 

3000.0 

410.0 

3800.0 

412.0 

4000.0 

400.0 

6500.0 

509.0 

6000.0 

472.0 

7800.0 

523.0 

7500  0 

457.0 

8000.0 

503.0 

8000.0 

475.0 

Profile  5 

Height  (ft) 

M-units 

0.0 

335.0 

400.0 

331.0 

2500.0 

406.0 

4500.0 

386.0 

7000.0 

476.0 

7500.0 

471  0 

8000.0 

489.0 

1  ransmitter  height  =  KOO  and  7000  It 
Antenna  beamwidth  -  0  5  and  15 
Antenna  elevation  .ingle  -  0  and  0  5 
Distance  between  profiles  -  50  mm 
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APPENDIX  A. 


Microsoft  QuickBASIC  (version  2.01)  program  code  for  use  on  IBM  compatible  PC, 
XT  or  AT  desktop  computer  with  MS-DOS  Version  3.2. 


A-i 


Program:  Laterally  Heterogeneous  Atmosphere  Raytrace 


Purpose:  Trace  a  ray  through  an  atmosphere  while  allowing  the 

refractive  index  to  vary  both  in  the  horizontal  and 
the  vertical. 

Origin:  Naval  Ocean  Systems  Center 

Code  542 

San  Diego,  CA  92152-5000 
Wayne  L.  Patterson 

Edition:  13  April  1987 

Glossary : 

a  -  array  index  counter 

a  -  intermediate  variable 

a$  -  temporary  string  variable 

angleinc  -  increment  of  beamwidth  (radians) 

b  -  intermediate  variable 

balpha  -  beginning  angle  (radians) 

beamwidth  -  antenna  beamwidth 

bht  -  beginning  height  (meters) 

brng  -  beginning  range  (kilometers) 

btmht  -  height  of  layer's  bottom  (meters) 

c  -  intermediate  variable 

ch  -  slope  of  the  layer  boundary 

charwd  -  number  of  screen  pixels  in  one  character  width 

charht  -  number  of  screen  pixels  in  one  character  height 

cm  -  slope  of  the  M-unit  gradient 

col  -  screen  column 

cursor  -  array  which  holds  the  image  of  a  plus  sign 

delh  -  change  of  height  (meters) 

delr  -  change  of  range  (kilometers) 

delu  -  change  of  M-units 

diff  -  difference  between  two  values 

dots  -  holder  of  refractive  condition  image 

dummy$  -  temporary  holder  for  operator  inputs 

e  -  vapor  pressure  (millibars) 

ealpha  -  ending  angle  (radians) 

eht  -  ending  height  (meters) 

elevangle  -  antenna  elevation  angle 

erng  -  ending  range  (kilometers) 

error?  -  conditional  checking  flag 

first  -  flag  to  signal  previous  use  of  program 

flag  -  temporary  counter 

fttom  -  feet  to  meters  conversion  factor 

h  -  temporary  height  (various  units) 

hiteO  -  height  at  level  0 

hitel  -  height  at  level  1 

hlim  -  height  of  tallest  profile  (meters) 

hmax  -  maximum  height  of  trace  (units  of  choice) 

ht  -  temporary  height 

htick  -  number  of  vertical  axis  tick  marks 
htop  -  temporary  height  (meters) 


hzero  -  radiosonde  launch  height 
i  -  array  index  counter 

j  -  array  index  counter 

k  -  array  index  counter 

layer  -  temporary  layer  counter 

lft  -  viewport's  left  screen  column 

lowest  -  lowest  height  within  a  profile  (meters  or 
pixels) 

lower  -  temporary  height  (meters) 

m  -  array  index  counter 

matkm  -  M-unit  value  at  height  of  one  kilometer 
matsfc  -  M-unit  value  at  the  surface 
mcolor  -  color  to  use  for  line  drawing 
mgrad  -  M-unit  gradient 

mitokm  -  miles  to  kilometers  conversion  factor 
mmax  -  maximum  M-unit  value  within  a  profile 

mmin  -  minimum  M-unit  value  within  a  profile 

mode$  -  flag  for  screen  type 

most  -  maximum  number  of  levels  within  all  profiles 
n  -  N-units 

n  -  array  index  counter 

ncol  -  temporary  screen  column  counter 

nmax  -  number  of  levels  within  each  environmental 

profile 

nnn  -  dimension  of  array  "dots" 

nprof  -  number  of  profiles 

nrays  -  number  of  rays  to  trace 

nrow  -  temporary  screen  row  counter 
offset  -  incremental  summation  of  rlim 
ok$  -  conditional  checking  flag 
op  -  option  number 

p  -  environmental  profiles 

index  1  -  counter  of  profile 
index  2  -  counter  of  levels 
index  3  -  environmental  data 

1  -  heights  (meters) 

2  =  M-units 

3  -  M-unit  gradient 

pix  -  same  a  p  but  data  given  in  screen  pixels 

pr  -  array  to  hold  profile  numbers  for  screen 

prompts 

pres  -  pressure  (millibars) 

presO  -  pressure  at  level  0 

presl  -  pressure  at  level  1 

putdots$  -  conditonal  checking  flag 

pzero  -  pressure  at  hzero 

rh  -  relative  humidity  (percent) 

rlim  -  range  between  profiles 

rngnow  -  temporary  range  (kilometers) 

rnginc  -  increment  of  range  (kilometers) 

row  -  screen  row 

startangle  -  ray's  initial  launch  angle  (radians) 
style%  -  style  to  use  for  line  drawing 
t  -  temporary  holder  of  environmental  data 


ta  -  temperature  (degrees  Kelvin) 

table  -  table  of  tick  mark  values  and  multiplication 
factors 

temp  -  temperature  (degrees  Celsius) 
title$  -  title  label 

tlevel  -  layer  containing  the  transmitter 
tmax  -  level  counter  for  array  t  (corresponds  to  nmax 
for  array  p) 

topht  -  height  of  layer's  top  (meters) 

tran  -  transmitter  height  (meters) 

tstrO  -  temperature  at  level  0 

tstrl  -  temperature  at  level  1 

units$  -  units  flag  for  English  or  metric 

upper  -  temporary  height  (meters) 

x  -  temporary  range  (km,  nmi,  or  pixels) 

xl  -  temporary  range  (km,  nmi.  or  pixels) 

xcol  -  screen  column  counter 

xinc  -  increment  of  range  (kilometers) 

xlabel$  -  horizontal  axis  labels 

xlim  -  maximum  range  of  trace  (kilometers) 

xlimit  -  temporary  range  limit  (kilometers) 

xmax  -  maximum  range  of  trace  (units  of  choice) 

xpixel  -  temporary  horizontal  screen  pixel 

xtick  -  number  of  horizontal  axis  tick  marks 

y  -  temporary  height  (m,  ft,  or  pixels) 

yl  -  temporary  height  (m,  ft,  or  pixels) 

ybtm  -  temporary  height  (pixels) 

yinc  -  increment  of  height  (meters) 

ylabel$  -  vertical  axis  labels 

ylimit  -  temporary  height  limit  (meters) 

ypixel  -  temporary  vertical  screen  pixel 

yrow  -  screen  row  counter 

ytop  -  temporary  height  (pixels) 


option  base  1 

dim  cursor#(26) ,  nmax(5),  p(5,150,3),  rlim(5) ,  pix(5 , 130 , 3) ,  s(5) 
dim  t(5, 150,2) 


common  shared  /environment/  runax(),  nprof,  p(),  units$, 

fttom,  mitokm,  t() 


common  shared  /tracecnst/  beamwidth,  elevangle,  hlim,_ 

nrays ,  rlim(),  tran 

common  shared  /graphcnst/  charwd,  charht,  cursor#(),  mode$,_ 

pix(),  x,  y,  hmax,  htick,_ 

xmax,  xtick,  s(),  putdots$ ,  xcol,  yrow 


common  shared  /errorhandel/  dummy$ 


nnn  -  15000 

dim  shared  dots (nnn) 


determine  what  type  screen  the  compute  is  using  and  set 
screen  type  and  character  size  accordingly 

charwd  -  8 
charht  -  14 


on  error  goto  cga:  screen  9 
on  error  goto  ega:  color  7 
mode?  -  "EGA16" 
goto  endset 

cga:  mode$  -  "CGA" 

screen  2 
charht  -  8 
resume  endset 

ega:  mode$  -  "EGA4" 

resume  endset 


endset:  on  error  goto  errorhandler 


print 
print 
print 
print 
print  " 
print  " 
print  " 
print  " 
print  " 
print 
print  ” 
print  " 
print  " 
print  " 


This  rayt race  program  will  allow  the  refractive" 
index  to  vary  horizontally  and  vertically.  Up  to" 
5  radiosonde  height  versus  M-unit  profiles  may  be" 
used  to  more  adequately  describe  the  horizontally" 
varying  atmosphere." 

For  further  information  contact  W.L.  Patterson," 
Naval  Ocean  Systems  Center,  Code  543,  San  Diego," 
CA,  92152-5000:  Commercial  619-225-7247:  Autovon" 
933-7247. " 


locate  25,1 

line  Input  "Press  <ENTER>  to  continue.  ";ok$ 
els 

ok$  -  "yes" 
first  -  0 
putdots$  -  "no" 

while  ok$  -  "yes" 
if  first  -  0  then 
call  envidata 
els 

call  mgradient 
call  drawprof 
call  sysdata 
call  drawrange 
call  raytrace 
first  -  1 
putdots$  -  "yes" 
else 

locate  25,1 

line  input  "Do  you  want  to  change  environmental  data?  (yes.no) 
ok$ 

if  lef t$ (ok$ , 1 )  -  "y"  or  left$(ok$,l)  -  "Y"  then 
putdots$  -  "no" 
call  envidata 
els 

call  mgradient 
call  drawprof 
call  sysdata 
call  drawrange 
call  raytrace 
putdots$  -  "yes" 

else 

call  sysdata 
call  drawrange 
call  raytrace 
end  if 
end  if 
locate  25,1 

line  input  "Do  you  want  to  rerun  the  program?  (yes, no)";  ok$ 
if  left$(ok$ , 1)  -  "y"  or  left$(ok$,l)  -  "Y"  then 
window 
view 
els 

ok$  -  "yes” 
end  if 
wend 
end 

'  subroutine  which  will  move  the  graphics  cursor  by  the  operator 

'  holding  the  arrow  keys 


downcursor : 

put  (x,y),  cursor#,  xor 
y  -  y  +  1 

if  y  >  13  *  charht  then  y  -  1 
put  (x,y),  cursor#,  xor 
return 
lef tcursor : 

put  (x,y),  cursor#,  xor 
x  -  x  -  5 


if  x  <  8  then  x  -  623 

put  (x,y) , 

cursor# , 

xor 

return 

rightcursor : 

put  (x,y), 
x  -  x  +  5 

cursor# , 

xor 

if  x  >  623 

then  x  - 

8 

put  (x , y) , 
return 

cursor# , 

xor 

upcursor : 

put  (x,y) , 

y  -  y  -  i 

cursor# , 

xor 

if  y  <0.5 

*  charht 

then 

put  (x,y) , 

cursor# , 

xor 

return 


*  charht 


'  error  handing  subroutines 

errorhandler : 

if  err  -  75  then  '  directory  all  ready  exists 

resume  next 

elseif  err  -  53  then  '  file  not  found 
close  #1 
beep 

locate  24,1 

print  "File  not  found.  Please  check  spelling  and  try  again, 
locate  25,1 

line  input  "Enter  the  desired  file  name.  ";dummy$ 
if  left$(dummy$ , 3)  -  "new"  or  lef t$ (dununy$ , 3)  -  "NEW"  then_ 
call  newdata 
resume  0 
end  if 

on  error  goto  0 
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Subrout in* 


envldata 


'  Purpose:  i.  Call  for  entry  of  environmental  data  from 

keyboard  or  read  the  data  from  a  previously 
created  disk  file. 

2.  Determine  the  x  and  y  axis  tick  mark  intervals 
based  upon  the  maximum  range  and  height  of 
the  data. 

'  3.  Insure  that  all  profiles  start  at  a  height  of 

zero  and  end  at  the  same  height 

sub  envidata  static 
dim  table(7,2) 

'  create  an  environmental  files  subdirectory  if  one  does 

'  not  all  ready  exist 


els  0 

mkdir  "enviros" 
start : 

print  "The  below  listed  environmental  files  may  be  used  for  data 
print  "input  by  typing  the  file  name  in  response  to  the  prompt.  If 
print  "new  data  is  desired,  respond  to  the  prompt  with  <new>." 
print 
print 

files  "enviros\" 
locate  25,1 

line  input  "Enter  the  desired  file  name.  dummy$ 
if  left$(dummy$ , 3)  -  "new”  or  lef t$ (dummy$ , 3)  -  "NEW"  then 
call  newdata 

else 

open  "enviros\"  +  dummy$  for  input  as  #1 
dummy$  -  input$ (42 7 , #1 ) 
units$  -  input$(3,#l) 
dummy$  -  input$(73 ,#1) 
input  #1,  nprof 
dummy$  -  input$(21 ,#1) 
f ttom  -  0 . 3048 
mitokm  -  1 . 854 
if  units$  -  "mks"  then 
fttom  -  1 
mitokm  -  1 
end  if 

for  i  -  1  to  5  '  read  number  of  levels 

input  #1,  nmax(i) 
next  i 

dummy$  -  input$ ( 39 , #1 ) 


for  i  -  1  to  4 

input  #1,  rlim(i) 


'  read  distance  between  profiles 


rlim(i)  -  rlim(i)  *  mitokm 
next  i 

dummy$  -  input$ ( 200 , #1 ) 
most  -  0 

for  i  -  1  to  nprof  '  determine  the  most  number  of  levels 

if  nmax(i)  >  most  then  most  -  nmax(i) 
next  i 

for  j  -  1  to  most  '  read  the  height  and  M-units 

for  i  -  1  to  nprof 

input  #1,  p ( i  ,  j  .  1 )  ,  p(i,j,2) 
p(i,j,l)  -  p(i,j,l)  *  fttom 
next  i 

dummv$  -  input$  ( 2  ,  #1 ) 
next  j 
close  #1 
end  if 


'  find  the  highest  height  of  all  the  profiles  and  the  total 

'  range 

hlim  -  0 
xlim  -  0 

for  i  -  1  to  nprof 

if  p( i , nmax( i) , 1)  >  hlim  then  hlim  -  p( i , nmax( i) , 1 ) 
if  i  <  nprof  then  xlim  -  xlim  +  rlim(i) 
next  i 

'  determine  the  number  and  intervals  of  the  vertical 

'  tick  mark 

restore  bound 
bound : 

data  4,4,5,5,8,4,10,5,15,5,20,4,25,5 
for  i  -  1  to  7 

read  table(i.l),  tabled, 2) 
next  i 

hmax  -  hlim  /  fttom 
factr  -  1 
aa : 

for  i  -  1  to  7 

if  hmax  <-  tabled, 1)  *  factr  +  0.1  then  goto  bb 
next  i 

factr  -  factr  *  10 
goto  aa 

bb: 

hmax  -  tabled, 1)  *  factr 
htick  -  tab  led, 2) 

determine  the  number  and  intervals  of  the  horizontal 


tick  marks 


xmax  -  xlim  /  mitokm 
factr  -  1 

cc : 

for  i  -  1  to  7 

if  xmax  <-  table(i.l)  *  factr  then  goto  dd 
next  i 

factr  -  factr  *  10 
goto  cc 

dd: 

xmax  -  table(i.l)  *  factr 
xtick  -  table(i,2) 

'  compute  the  M-unit  value  at  one  kilometer  for  each  profile 

for  i  -  1  to  nprof 
ok$  -  "no" 
j  -  nmax(i) 

if  p(i,j,l)  <  1000  then 

diff  -  (1000  -  p(i,j ,1))  /  1000 
matkm  -  (118  *  diff)  +  p(i,j,2) 
elseif  p(i,j,l)  -  1000  then 
matkm  -  p( i , j , 2 ) 

else 

while  ok$  -  "no" 

if  p(i,j,l)  >  1000  then 
J  -  j  -  1 

else 

matkm  -  p(i,j+l,2)  -  (p(i,j+l,2)  -  p ( i , j , 2 ) )  *_ 
(p(i,j+l,l)  -  1000)  /  _ 

(p(i,j+l.l)  -  p(i.j.l)) 
ok$  -  "yes" 
end  if 

wend 
end  if 

if  necessary,  compute  the  M-unit  value  at  the  height  (hlim) 
'  and  insert  it  into  the  profile 

if  p( i , nmax( i ) , 1)  <  hlim  then 
nmax(i)  -  nmax(i)  +  1 
p ( i , nmax( i ) , 1 )  -  hlim 
h  -  hlim  /  1000 

p(i ,nmax( i) ,2)  -  157  *  h  +  p(i,l,2)  *  exp ( - log(p ( i , 1 , 2 )_ 

/  (matkm  -  157))  *  h) 

end  i  f 
next  i 
end  sub 


A9 


O  K.^  O  A.", 


Subroutine:  sysdata 


Purpose:  Solicit  the  transmitter  data  from  the  operator 


sub  sysdata  static 
els  0 

ok$  -  "no" 
while  ok$  -  "no" 

if  units$  -  "fps"  then 
locate  24 , 1 

print  "Transmitter  must  be  between  0  and";_ 
int(hlim/fttom) ;  "feet.” 
locate  25,1 

line  input  "Enter  transmitter  height  (feet)  ";  dummy$ 
else 

locate  24,1 

print  "Transmitter  must  be  between  0  and";_ 
int(hlim) ;  "meters." 
locate  25,1 

line  input  "Enter  transmitter  height  (meters)  ";  dummy$ 
end  if 

tran  -  val(dummy$)  *  fttom 
if  tran  <  0  or  tran  >  hlim  then 
beep 

print  "Transmitter  height  must  be  between  0  and 
hlim  /  fttom;".  Try  again." 

else 

ok$  -  "yes" 
end  if 
wend 

els  0 

locate  25,1 
ok$  -  "no" 
while  ok$  -  "no" 

input  "Enter  beamwidth  in  degrees  (0.5  to  45).  ",  beamwidth 
if  beamwidth  <0.5  or  beamwidth  >  45  then 
beep 

print  "Beamwidth  must  be  between  0.5  and  45.  Try  again." 
else 

ok$  -  "yes" 
end  if 
wend 

els  0 

locate  25,1 
ok$  -  "no" 
while  ok$  -  "no" 

input  "Enter  elevation  angle  In  degrees  (-10  to  10).  ", 
elevangle 

if  elevangle  <  -10  or  elevangle  >  10  then 
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if 
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print  "Elevation  angle  must  be  between  -10  and  10.  Try  again." 
else 

ok$  -  "yes" 
end  if 
wend 

els  0 

locate  25,  1 
ok$  -  "no 
while  ok$  -  "no" 

input  "Enter  number  of  rays  to  trace  (1  to  100).  ",  nrays 
if  nrays  <  1  or  nrays  >  100  then 
beep 

print  "Number  of  rays  must  be  between  1  and  100.  Try  again." 
else 

ok$  -  "yes" 
end  if 
wend 
els  0 
end  sub 
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Subroutine:  drawprof 


'  Purpose:  1.  Draw  the  M-unit  versus  height  profiles  on  the 

'  screen 

'  2.  Determine  if  layers  can  be  automaticly  defined 

'  If  not,  call  for  layer  drawing  subroutine. 

sub  drawprof  static 

'  print  the  graph  title  and  vertical  axis  label 

locate  1,21 

print  "SELECTED  HEIGHT  VERSUS  M  UNIT  PROFILES" 
a$  -  "HEIGHT" 
row  -  9 

for  i  -  1  to  6 
locate  row,l 
print  mid$(a$,i,l) 
row  -  row  +  1 
next  i 

'  define  the  left  margin  (column  number)  of  each  drawing  cell 

restore  bounds 

bounds:  data  3,  19,  35,  51,  67 
col  -  6 

'  loop  to  draw  each  profile  cell 

for  i  -  1  to  nprof 

'  print  cell  title  and  horizontal  axis  label 

locate  3,col-l 

print  "PROFILE"  +  str$(i) 

locate  20, col 

print  "M  UNITS" 

col  -  col  +  16 

open  the  viewport  and  find  the  mininum  and  maximum 
'  x  axis  value 

read  1ft 

view  (lft*charwd,4*charht)  -  ((lft+ll)*charwd,  18*charht) 
mmin  -  1 . e5 
mmax  -  0 

for  j  -  1  to  nmax(i) 

if  p ( i , j , 2 )  >  mmax  then  mmax  =  p ( i , j , 2 ) 

if  p ( i , j , 2 )  <  mmin  then  mmin  =  p(i,j,2) 

next  j 

'  define  the  drawing  window  and  draw  the  horizontal  axis 

window  (mmin,  p ( i , nraax ( i ) , 1 )/  fttom)  -  (mmax,  0) 
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line  (mmin, 0)  -  (mmax,0),14 

draw  the  m  unit  versus  height  profile  using  various  colors 
for  the  gradients 

for  j  -  1  to  nmax(i)  -  1 
if  p(i, j ,3)  <=  0  then 
mcolor  -  4 
style%  -  &hffff 
elseif  p(i , j ,3)  <  7.9e-5  then 
mcolor  -  5 
style%  -  &hfl8f 

elseif  p( i , j ,3)  <-  1.57e-4  then 
mcolor  -  2 
style%  -  &hcccc 

else 

mcolor  -  1 
style%  -  &hff00 
end  if 

if  mode$  O  "CGA"  then  style%  -  &hffff 
line  (p(i,j,2),  p ( i , j , l)/fttom) 

(p(i.j+1.2),p(i,j+l, 1) /ft tom) 
mcolor,  ,  style% 

fill  the  pix  array  with  x  axis  pixels,  y  axis  heights 
and  z  axis  l.e7 

pix(i,j+l,l)  -  point(0)+lft*charwd 
pix(i,j+l,2)  -  point(l) 

if  p(i,j+l,l)  -  0  then  pix(i,j+l,3)  -  l.e7 
if  p(i,j+l,l)  O  0  then  pix(i,j+l,3)  -  p(i,j+l,l) 
next  j 

next  i 

window 
view 
row  -  22 

call  refcond(row) 

check  to  see  if  all  profiles  have  same  number  of  points 

ok$  —  "yes" 
i  -  1 

while  ok$  -  "yes"  and  i  <«  nprof 
if  nmax(l)  o  nmax(i)  then 
ok$  -  "no" 
flag  -  -1 

else 

i  -  i  +  1 
end  if 

wend 


if  all  profiles  have  same  number  of  points,  count  number  of 


levels  on  each  profile  that  has  a  negative  m-unit  gradient 


9-J 

D 


..‘•Vtl 


a  a 

!jl 

ij 

lyjj 

I 


if  ok$  -  "yes"  then 
i  -  1 

while  i  <—  nprof 
flag  -  0 

j  “  1 

while  j  <=  nmax(l) 

if  p(i,j,3)  <  0  then  flag  -  flag  +  1 

j  -  j  +  1 

wend 

i  -  i  +  1 

wend 
end  if 

’  draw  boundary  lines  and  ask  for  approval  if  equal  number  of 

'  points  otherwise  go  straight  to  operator  drawing  routine 

if  ok$  -  "yes"  then 
window 

view  (1,  4*charht)  -  (639,  18*charht) 

mcolor  -  4 

for  j  -  2  to  nmax(l) 

for  i  -  1  to  nprof  -  1 

if  mode$  -  "EGA16"  then 

if  p(i,j,3)  <-  0  then 
mcolor  -  4 

elseif  p(i , j ,3)  <  7.9e-5  then 
mcolor  -  5 

elseif  p(i, j ,3)  <-  1.57e-4  then 
mcolor  -  2 

else 

mcolor  -  1 
end  if 
end  if 

line  (pix(i, j ,1) ,pix(i, j ,2))  -  (pix(i+l , j , 1) 
pix( i+1 , j , 2) ) ,  mcolor,  ,  &hffff 

next  i 
next  j 
locate  23,1 
beep 

if  flag  -  0  then 

print  "Note!  Trapping  layer  not  continuous  between  profiles." 

else 

print  "Note!  Multiple  trapping  layers  within  a  profile." 
end  if 

line  input;  "Want  to  define  the  layers  yourself?  (yes  or  no) 
ok$ 

if  left$(ok$ , 1)  -  "y"  or  left$(ok$,l)  -  "Y"  then 

erase  the  boundary  line  previously  drawn  and  call  for  drawing 
'  routine 
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for  j  =  2  to  nmax(l) 

for  i  =  1  to  nprof  -  1 

line  (pix(i,j,l),pix(i.j,2))  -  (pix(i+l,j,l) 
pix(i+l,j,2)),  0,  ,  &hffff 


next  i 


next  j 
window 


call  drawlayers 
end  if 


window 


beep 
if  flag 


if  flag  =  -  1  then 
locate  23,1 

print  "Profiles  contains  unequal  number  of  levels." 
end  i f 
locate  24 , 1 

line  input;  "Press  <ENTF.R>  to  continue.  ok$ 
call  drawlayers 
end  if 


els  0 
end  sub 


Subroutine:  drawlayers 


'  Purpose:  1.  Prompt  the  operator  for  determination  of  layer 

'  boundaries. 

'  2.  Re-establish  the  profiles  by  including  any  new 

'  layers  as  defined  by  the  operator 

sub  drawlayers  static 

'  fill  the  prompt  array.  These  prompts  are  profile  numbers  which 

'  will  be  considered  next. 

dim  pr(5,5,5),  tmax(150) 

restore  a 

for  i  -  1  to  5 

for  j  -  1  to  5 

for  k  -  1  to  5 

read  pr(i , j ,k) 
next  k 
next  j 
next  i 


data  0,0, 0,0, 0.0, 0,0, 0,0, 0,0, 0,0, 0,0,0, 0,0, 0,0, 0,0, 0,0 
data  2, 0,0, 0,0, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 
data  2, 3, 0,0, 0,3, 2, 1,0, 0,2, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 
data  2, 3, 4, 0,0, 3, 4, 3, 1,0, 4 , 2, 2, 1,0, 3, 2, 1,0, 0,0, 0,0, 0,0 
data  2, 3, 4, 5, 0,3, 4, 5, 4, 1,4, 5,3, 2, 1,5, 2, 3, 2, 1,4, 3, 2, 1,0 

'  fill  the  first  layer  of  the  temporary  profile  array 

'  with  the  first  layer  of  the  environmental  array 

for  i  -  1  to  nprof 

t(i . 1 , 1)  -  p( i , 1 ,  1 ) 
t ( i  ,  1 , 2 )  =  p( i , 1 , 2) 
tmax(i)  =  1 
next  i 


'  define  the  graphics  crosshair  and  save  the  image 

view(l,  21  *  charht)  -  (639,  23  *  charht) 

els 

view 

locate  2  , 2 

line  (1,  0.5*  charht)  -  (2  *  charwd,0.5  *  charht),  14 
line  (charwd,  1)  -  (charwd,  charht),  14 
get  (1,  1)  -  (2  *  charwd,  charht),  cursor# 
put  (1,  1),  cursor#,  xor 

view  (1,  4  *  charht)  -  (639,  18  *  charht) 

find  the  most  number  of  points  in  any  profile 
most  -  0 
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to i  i  -  1  to  nprof 

If  nmax(i)  >  most  then  most  -  nmax(i) 


'  start  the  loop  for  solicitation  of  points 

c:  k  $  -  "yes" 

layer  -  2 

while  ok$  O  "no" 

find  a  starting  point  -  the  lowest  nonused  point  of  all  profiles 
ok$  -  "no" 

j  ~  2 

lowest  -  0 

while  ok$  -  "no"  and  j  <-  most 
for  i  -  1  to  nprof 

if  pix(i,j,3)  O  l.e7  and  pix(i,j,2)  >  lowest  then 
n  -  i 

jj  "  j 

lowest  -  pix(i , j , 2) 
ok$  -  "yes" 
end  if 
next  i 

if  ok$  -  "no"  then  j  -  j  +  1 

wend 

if  ok$  -  "yes"  then 
k  -  1 

for  j  -  1  to  nprof  -1 

put  the  starting  point  into  the  temporary  profile  array 
and  flag  it  as  used 

t(n, layer, 1)  -  p(n,jj,l) 
t(n,layer,2)  -  p(n,jj,2) 
pix(n, j j ,3)  -  1 . e7 

put  the  crosshair  at  the  starting  point 

i f  j  -  1  or  k  -  pr (nprof , n , k)  then 

the  next  if  statement  is  a  special  case  fix 

if  nprof  -  4  and  n  -  4  and  k  -  2  then  goto  10 
x  -  pix(n,jj,l)  -  charwd 
y  -  pix(n,jj,2)  0.5  *  charht 
xl  -  pix(n,jj,l) 
yl  -  pix(n, j j ,2) 
put  (x,y),  cursor**,  xor 
if  )  O  1  then  k  -  k  +  1 
end  i  f 
lU 

m  •*  pr  (nprof  ,  n  ,k) 


error$  -  "yes" 
while  error?  -  "yes 


i 

$ 
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'  solicit  for  next  point 

locate  24,1 
beep 

print  "Using  arrow  keys,  select  point  from  profile" 
pr(nprof ,n,k) ;  "then  press  ENTER. 

'  turn  on  the  arrow  keys  and  event  trap 

key  (11)  on 
key  (12)  on 
key  (13)  on 
key  (14)  on 

on  key  (11)  gosub  upcursor 
on  key  (12)  gosub  leftcursor 
on  key  (13)  gosub  rightcursor 
on  key  (14)  gosub  downcursor 

incursor : 

a$  —  " " 
while  a$  -  "" 
a$  -  inkey$ 

wend 

if  a$  -  chr$(13)  then  goto  nextprof 
goto  incursor 

nextprof : 

'  convert  the  input  y  pixel  to  height 

h  =  y  +  0. 5  *  charht 

ytop  -  1 

ybtm  =  14  *  charht 
htop  -  p(m,nmax(m) , 1) 

h  =-  htop  -  htop  *  (ytop-h)  /  (ytop-ybtm) 

'  find  the  lowest  nonused  point  (a)  on  the  next  profile 

a  -  2 

for  i  -  1  to  nraax(m) 

if  pix(m,i,3)  -  l.e7  then  a  -  a  +  1 
next  i 

’  compute  acceptable  bounds  about  the  point  "a" 

upper  -  p(m,a,l)  +  p(m,a,l)  *  0.03 
lower  -  p(m, a , 1)  -  p(m,a,l)  *  o.05 

'  check  the  input  height's  relationship  to  the  bounds 


if  h  >  t(m, layer- 1 , 1)  then 


input  height  less  then  bottom  bound  and  greater  than  the  last 
temporary  profile  point  so  interpolate  for  M-unit  value  and 
insert  the  new  level  into  the  temporary  profile  array 

t(m, layer , 2)  -p(m,a,2)  -  ((p(m,a,2) 

p(m,a-l,2))  *  (p(m, a , 1)  -  h)_ 

/  (p(m, a, 1)  -  p(m, a-1 , 1) ) ) 
t (m, layer , 1)  -  h 
tmax(m)  -  tmax(m)  +  1 
tmax(n)  -  tmax(m) 
error$  -  "no" 

else 

input  height  less  than  previous  temporary  profile  point  so  give 
crossing  bounds  error  message 

beep 

locate  23,1 

print  "Boundaries  can't  cross.  Press  ENTER."; 


pausea : 

a$  - 

while  a$  -  "" 
a$  -  inkey$ 

wend 

if  a$  O  chr$(13)  then  goto  pausea 
locate  23,1 
print 
end  if 

elseif  h  >  upper  then 

'  input  height  greater  than  bounds  so  give  missed  breakpoint  error 

'  message 


beep 

locate  23,1 

print  "M-unit  breakpoint  exceeded.  Press  ENTER." 


pauseb : 

a$  -  "" 
while  a$  -  "" 
a$  -  inkey$ 

wend 

if  a$  O  chr$(13)  then  goto  pauseb 

locate  23,1 

print 

else 

input  height  within  bounds  so  insert  the  point  (a)  into  the 
'  temporary  profile  array  and  flag  it  as  being  used 


t(m, layer , 2)  -  p(m,a,2) 
tmax(m)  -  tmax(m)  +  1 
tmax(n)  -  tmax(m) 
pix(m,a,3)  -  l.e7 
error$  -=  "no" 
end  if 

if  error$  -  "yes”  then 

if  there  was  an  error,  move  crosshairs  back  to  starting  point 

put  (x,y),  cursor#,  xor 
x  -  xl  -  charwd 
y  -  yl  -0.5  *  charht 
put  (x,y),  cursor#,  xor 


else 


Eft 


if  no  error,  draw  boundary  line  between  the  two  profiles 
erase  the  crosshairs  if  on  the  last  profile  and  move  on 
to  the  next  profile 


line  (xl,  yl)  -  (x  +  charwd,  y  +0.5  *  charht),  14 
xl  -  x  +  charwd 
yl  -  y  +  0.5  *  charht 
if  j  -  nprof  -1  then 

put  (x,y),  cursor#,  xor 
elseif  j  -  nprof  -  n  then 
put  (x,y),  cursor#,  xor 
end  if 
k  -  k  +  1 
end  if 

wend 
next  j 

layer  -  layer  +  1 
end  if 


wend 


refill  the  profile  array  with  the  temporary  height  M-unit  values 
and  recompute  the  M-unit  gradients 


for  i  -  1  to  nprof 

tmax(i)  -  tmax(i)  +  1 
t( i , tmax( i ) , 1 )  -  p( i , nmax( i) , 1 ) 
t(i , tmax( i) , 2)  -  p( i , nmax( i) , 2) 
next  i 


for  i  -  1  to  nprof 

nmax(i)  -  tmax(i) 
for  j  -  1  to  nmax(i) 

p(i,j  ,1)  -  t ( i , j ,1) 
P ( i . j  ,2)  -  t( i , j ,2) 
next  j 
next  i 

call  mgradient 
end  sub 


!w!w! 


Subroutine:  Drawrange 


*  Purpose:  1.  Draw  the  axis  and  labels  for  the  raytrace 

sub  drawrange  static 

xlabel$  -  "Range  in  Nautical  Miles" 
ylabel$  -  "Height  in  Feet" 

title$  -  "Raytrace  -  Laterally  Heterogeneous  Atmosphere" 
if  units$  -  "mks"  then 

xlabel$  -  "Range  in  Kilometers" 
ylabel$  -  "Height  in  Meters" 
end  if 


'  establish  drawing  limits  and  tick  mark  increments 

xlimit  -  xmax 
yllmit  -  hmax 
xinc  -  xlimit  /  xtick 
yinc  -  ylimit  /  htick 

'  determine  the  start  points  based  upon  the  terminal  screen  being 

'  used  and  draw  the  axis 

view 


if  htick 

- 

A  then 

yrow 

- 

18.5 

nrow 

- 

A 

else 

yrow 

- 

17.5 

nrow 

- 

3 

end  if 

if  xtick 

- 

A  then 

xcol 

- 

76.5 

ncol 

- 

16 

else 

xcol 

- 

72.5 

ncol 

- 

12 

end  if 

line  (12.5  *  charwd,  2.5  *  charht)  -  (12.5  *  charwd,  yrow  *  charht),7 
line  (12.5  *  charwd,  yrow  *  charht)  -  (xcol  *  charwd,  yrow  *  charht), 7 
ypixel  -  2.5  *  charht:  xpixel  -  10  *  charwd 

'  draw  and  label  the  vertical  axis  tick  marks 

row  -  3 
fac  -  htick 

for  i  -  2  *  htick  to  0  step  -1 
if  i  mod  2-0  then 

line  (xpixel , ypixel)  -  (xpixel  +  2.5  *  charwd,  ypixel) 
locate  row,  A 
print  yinc  *  fac 
fac  -  fac  -  1 
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row  -  row  t-  nrow 


else 

line  (xpixel+charwd ,  ypixel)  -  (xpixel+2 . 5*charwd ,  ypixel) 
end  if 

ypixel  -  ypixel  +  nrow  *  charht  /  2 
next  i 

'  draw  and  label  the  horizontal  axis  tick  marks 

col  -  12 
fac  -  0 

ypixel  -  yrow  *  charht 
xpixel  -  xpixel  +  2.5  *  charwd 
for  i  -  0  to  2  *  xtick 
if  i  mod  2  “  0  then 

line  (xpixel,  ypixel)  -  (xpixel,  ypixel  +  charht) 

locate  20,  col 

print  xinc  *  fac 

fac  -  fac  +  1 

col  -  col  +  ncol 

else 

line  (xpixel,  ypixel)  -  (xpixel,  ypixel  +0.5  *  charht) 
end  if 

xpixel  -  xpixel  +  ncol  *  charwd  /  2 
next  i 

'  print  the  title  and  horizontal  axis  labels 

locate  1,  (80  -  len(title$))  /  2:  print  title? 
locate  21,  (80  -  len(xlabel$) )  /  2:  print  xlabel$ 

'  print  the  vertical  axis  label 


for  i  -  1  to  len(ylabel$) 

locate  (20  -  len(ylabel$ ) )  /  2  +  i,  2 
print  mid$ (ylabel$ , i , 1) 
next  i 

if  screen  is  enhanced  graphics  with  extra  memory,  call  for 
'  refractive  conditions  legion 

row  -  22 

if  mode$  -  "EGA16"  then  call  refcond(row) 

'  define  the  view  port  and  window  for  the  raytrace  graph 

view  (12.5  *  charwd,  2.4  *  charht)  -  (xcol  *  charwd,  yrow  *  charht) 
window  (O.hmax)  -  (xmax.O) 
end  sub 
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Subroutine:  refcond 

Purpose:  1.  Draw  the  refractive  condition  legion  at  the 

bottom  of  the  screen 


sub  refcond(row)  static 

locate  row, 10 

print  "TRAPPING" 

locate  row, 27 

print  "SUPERREFRACTIVE" 

locate  row, 51 

print  "STANDARD" 

locate  row, 68 

print  "SUBREFRACT I VE" 

'  define  the  box  boundaries,  interior  colors  and  line  style 

restore  limits 

limits:  data  3,  4,  &hffff,  20,  5,  &hfl8f,  44,  2,  &hcccc ,  61,  1,  &hff00 

'  read  the  drawing  attributes,  draw  and  fill  the  box 

for  i  -  1  to  4 

read  1ft,  clr,  style% 
if  mode$  -  "EGA16"  then 

view  screen  (lft*charwd,  (row-l)*charht)  - 
( (lft+5)*charwd,  (row)*charht) 
paint  ( (lft+l)*charwd,  (row-0. 5) *charht) ,  clr,  clr 

else 

line  (lft*charwd,  row*charht  -0.5*charht)  - 

( (lft+5)*charwd,  row*charht  -0. 5*charht) ,  clr,  ,  style% 


end  if 


Subroutine:  newdata 


'  Purpose:  1.  Solicit  the  operator  for  environemntal  data. 

'  2.  If  data  entered  in  units  other  than  M-units, 

'  convert  to  M-units 

sub  newdata  static 

dim  e(31),  pres(31),  rh(31),  table(7,2),  temp(31) 

'  initialize  arrays  and  constants 

erase  e  nmax  p  pres  rh  rlim  temp 

units$  -"fps" 

f ttora  -  0 . 3048 

mitokm  -  1.854 

els 

locate  4,1 
print 

print  "Atmosphere  specification  options  are:” 
print 

print  "  1.  Pressure,  temperature,  and  relative  humidity" 

print  "  2.  Height  and  N-units" 

print  "  3.  Height  and  M-units" 

ok$  -  "no" 
while  ok$  -  "no” 
locate  25,1 

input;  "Enter  specification  option  (1,  2  or  3)?  ",  op 
if  op  <  1  or  op  >  3  then 
beep 

print  "Option  must  be  between  1  and  3.  Try  again." 

else 

ok$  -  "yes" 
end  if 

wend 

els 

locate  25,1 

line  input;  "Enter  units  of  height/range  (english  or  metric). 
dummy $ 

if  left$(dummy$ , 1)  -  "m"  or  lef t$ (dummy$ , 1)  -  "M"  then 
fttom  -  1 
mitokm  -  1 
units$  -  "mks" 
end  if 

els 

ok$  -  "no" 
while  ok$  -  "no" 
locate  25,1 

input  "How  many  profiles  do  you  want  to  enter  (1  to  5)?  ",  nprof 
if  nprof  <  1  or  nprof  >  5  then 
beep 


I 


print  "Number  of  profiles  must  be  between  1  and  5.  Try  again. 


else 

ok$  -  "yes" 
end  if 


wend 


enter  the  atmosphere  data 


for  i  -  1  to  nprof 
els 

locate  1 , 1 

print  "Data  entry  for  profile";  i;  "of";  nprof 
print 

J  -  1 

ok$  -  "no" 

while  ok$  -  "no"  and  j  <-30 
print  "For  level";  j 
if  op  -  1  then 

if  j  -  1  then 

if  units$  -  "fps"  then 

input  "Enter  radiosonde  launch  ht  (ft),  pres  (mb)_ 
temp  (C)  and  rh(%).  ",  hzero,  pres(j),  temp(j),  rh(j) 
else 

input  "Enter  radiosonde  launch  ht  (meters) ,  pres_ 
(mb),  temp(C)  and  rh(%) .  ".hzero,  pres(j),  temp(j),  rh(j) 
end  if 

p ( i  ,  j , 1 )  -  hzero  *  fttom 
nmax ( i )  -  nmax ( i )  +1 
j  -  j  +  1 

else 

input  "Enter  pressure  (mb),  temp  (C) ,  and  rh  (-1,  -1,_ 
-1  to  end)  ",  pres(j),  temp(j),  rh(j) 
if  pres(j)  <  0  then 
ok$  -  "yes" 

elseif  pres(j)  >-  pres(j-l)  then 
beep 

print  "Pressure  must  decrease.  Try  again." 

else 

nraax(i)  -  nmax(i)  +  1 

j  -  j  +  1 

end  if 
end  if 

else 

if  units$  -  "mks"  then 
if  op  -  2  then 

input  "Enter  height  (meters)  and  N  units  (-1,-1  to 

ht ,  n 

m  -  n  +  ht/6.371 

else 

input  "Enter  height  (meters)  and  M  units  (-1,-1  to 

ht ,  m 

end  i  f 

else 

if  op  -  2  then 
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end)  ",  ht ,  n 


input  "Enter  height  (feet)  and  N-units  (-1,  -1  t 
m  -  n  +  ht*fttom/6 . 371 


end) 


else 

input  "Enter  height  (feet)  and  M-units  (-1.-1  to 

ht ,  m 

end  if 
end  if 

ht  -  ht  *  fttora 
if  ht  >—  0  then 
if  j  -  1  then 

P(i  .  j .1)  “  ht 
p(i,  j  ,2)  -  m 
nmax ( i )  -  j 

j  -  j  +  1 

else 

if  ht  <-  p(i,j-l,l)  then 
beep 

print  "Heights  must  be  in  increasing  order. 


Try  again.' 


else 

P(i.j.l)  “ 
P(i, j ,2)  - 
nmax ( i )  -  j 

j  -  j  +  1 

end  if 
end  if 


ht 

m 


else 

ok$ 

end  if 
end  if 


'yes' 


wend 


if  pressure  entered,  convert  to  height  and  compute  M-units 

if  op  -  1  then 

hiteO  -  hzero 
presO  -  pres(l) 
tstrO  -  temp(l)  +  273.2 
for  j  -  1  to  nmax(i) 

ta  -  temp(j)  +  273.2 
presl  -  pres(j) 

ee  *  6.105  *  exp(25.22  *  ( ta-273 . 2)/ta  -  5.31  *  _ 
log(ta/273 . 2) )  *  rh(j)/100 
tstrl  -  ta  +  .3794017  *  ta  *  ee/  (presl  -  ee) 
hitel  -  hiteO  +  14.643  *  (tstrl  +  tstrO)  * 
log(presO/presl) 
p(i , j , 1)  -  hitel 

p ( i , j , 2 )  -  77,6/ta  *  (presl  +  4810  *  ee/ta)  +_ 
hitel/6 .371 


end  if 


'  if  necessary,  set  surface  height  to  zero  and  compute  a  M-unit 

'  value 

if  p(i,l,l)  O  0  then 

for  j  -  nmax(i)  +  1  to  2  step  -1 
p(i,j ,1)  -  p(i,J-l.l) 
p(i,j ,2)  -  p(i, j -1.2) 
next  j 

P(i.l.l)  -  0 

p(  i  ,  1 , 2 )  -  p(  i ,  1 , 2)  -  lj.8  *  p(i,2,l)/1000 
nmax(i)  -  nmax(i)  +  1 
end  if 

els  * 

locate  25,1 
if  nprof  -  1  then 

if  units$  **  "fps"  then 

input  "Enter  range  (nautical  miles)  for  raytrace. 
rlim( i) 

else 

input  "Enter  range  (kilometers)  for  raytrace.  ",  rlim(i) 
end  if 

rlim(i)  -  rlim(i)  *  mitokm 
for  j  -  1  to  nmax(i) 

p(2 , j ,1)  -  p(l, j  ,1) 

P ( 2  ,  j ,2)  -  p(l , j ,2) 
next  j 

nmax ( 2 )  -  nmax ( i ) 
nprof  -  2 

elseif  i  <  nprof  then 

if  units$  -  "fps"  then 

input  "Enter  range  (nautical  miles)  to  next  profile. 
rlim( i) 

else 

input  "Enter  range  (kilometers)  to  next  profile.  ",_ 
rlim(i) 

end  if 

rlim(i)  -  rlim(i)  *  mitokm 
end  if 
next  i 

'  solicit  the  operator  for  a  storage  file  name  if  desired 

els 

ok$  -  "no" 
while  ok$  -  "no" 
locate  25,1 

line  input;  "Do  you  want  to  store  this  data  for  future  use  _ 
(yes.no)?",  dummy $ 

if  left$(duraray$ , 1)  -  "y"  or  left$ (dummy$ , 1)  -  "Y"  then 
els 
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locate  25,1 

line  input;  "Enter  file  name  for  storage  of  data  (8  characters_ 
max) .  " ;  dummy? 

bpen  "enviros\"  +  dummy?  for  output  as  #1 

print  #1,  "This  file  contains  environmental  data  for  the  multi 
-profile" 

print  #1,  "raytrace.  If  you  edit  this  file,  the  colons  (:)  mus_ 
t  remain" 

print  #1,  "in  place.  All  numbers  must  be  separated  by  at  1 
east  one" 

print  #1,  "space.  The  column  structure  of  the  heights  and  M-un 
its  must" 

print  #1,  "be  retained.  Heights  must  be  in  increasing  order. 
There  is" 

print  #1  *  "no  error  checking  on  a  stored  file. 

« 

print  #1,  " 

It 

print  #1,  using  "\  \" ;  units?; 

print  #1,  using  ; "  :units  of  height/range . _ 

mks  -  metric.  fps  =-  English." 

print  #1,  using  "#" ;  nprof; 

print  #1,  using  ; "  : number  of  profiles  " 

most  -  0 
for  i  -  1  to  5 

print  #1,  using  "###";  nmax(i); 
if  nmax(i)  >  most  then  most  -  nmax(i) 
next  i 

print  #1,  using  "  : number  of  levels  within  each  profi 

le" 

for  i  -  1  to  4 

print  #1,  using  "#####";  rlim(i)  /  mitokm; 
next  i 

print  #1,  using  .'distance  between  profiles  in  units  of  a 

bove" 

print  #1,  "" 

print  #1,  ”  Profile  1  Profile  2  Profile  3  Prof 

le  4  Profile  5” 

print  #1,  "  height  M-unit  height  M-unit  height  M-unit  height 
t  M-unit  height  M-unit" 
print  #1,  "" 

for  j  -  1  to  most 

for  i  -  1  to  nprof 

print  #1,  using  "#####.#";  p(i,j,l)  /  fttom; 
print  #1,  using  "  ####.#  ";  p ( i , j  , 2 ) ; 
next  i 

print  #1,  chr?(13) 
next  j 
close  #1 


Subroutine:  raytrace 


'  Purpose:  1.  By  taking  a  range  step,  compute  an  ending 

'  height  and  angle  for  a  ray  with  a  specified 

'  launch  angle . 

sub  raytrace  static 

dim  offset(5) ,  ch(5,150),  cm(5,150) 

startangle  -  (elevangle  -  (0.5  *  beamwidth))  *  0.0174532 
angleinc  -  (beamwidth  /  nrays)  *  0.0174532 

'  total  the  range  increments  and  compute  incremental  range  offsets 

xlimit  -  0 

for  i  »  1  to  nprof  -  1 

xlimit  —  xlimit  +  rlim(i) 
next  i 

offset(l)  -  0 

offset(2)  —  rlim(l) 

offset(3)  -  offset(2)  +  rlim(2) 

offset(4)  -  offset(3)  +  rlim(3) 

offset(5)  -  offset(4)  +  rlim(4) 

compute  the  layer  height  and  M-unit  gradient  coefficients 

for  i  -  1  to  nprof  -  1 
for  j  -  1  to  nmax(i) 

ch(i,j)  -  (p(i+l, j , 1)  -  p(i , j , 1) )  /  rlim(i) 
cm(i,j)  -  (p(i+l , j , 3)  -  p(i , j , 3) )  /  rlim(i) 
next  j 
next  i 

'  ask  for  display  of  refractive  conditons 

if  putdots$  -  "no"  then 
locate  25,1 

if  mode$  -  "EGA16"  then 

line  input;  "Want  to  see  all  conditions  or  just  trapping  (all 
or  trap)?  ",dummy$ 

else 

print  " 

tl  . 

end  if 

'  draw  dots  for  refractive  conditions 

locate  25,1 

if  mode$  O  "EGA16"  then 

print  "Drawing  trapping  regions.  Please  standby. 


print  " 

H  . 

end  if 

if  xtick  -  4  then  xinc  -  xlimit  /  128 
if  xtick  -  5  then  xinc  -  xlimit  /  120 
if  mode$  O  "CGA"  then 

if  htick  -  4  then  yinc  -  int(hlim)  /  112 

if  htick  -  5  then  yinc  -  int(hlim)  /  105 

else 

if  htick  -  4  then  yinc  -  int(hlim)  /  64 

if  htick  -  5  then  yinc  -  int(hlim)  /  60 

end  if 
y  -  yinc 
i  -  1 

j  -  1 

n  -  1 

while  y  <  hlim  -  yinc 
x  -  xinc 

while  x  <  xlimit 

ytop  -  ch(i , j+1)  *  (x-offset(i) )  +  p(i,j+l,l) 
ybtm  -  ch(i, j )  *  (x-offset(i) )  +  p(i,j,l) 
if  y  >  ytop  then 

j  -  J  +  1 

elseif  y  <  ybtm  then 

j  -  j  -  1 

end  if 

mgrad  -  cm(i,j)  *  (x-offset(i) )  +  p(i,j,3) 
if  mode$  -  "EGA16"  and  ( left$ (dummy$ , 1 )  -  "a"  or 
left$(dummy$ , 1)  -  "A")  then 
if  mgrad  <-  0  then 
mcolor  -  4 

elseif  mgrad  <  7.9e-5  then 
mcolor  -  5 

elseif  mgrad  <-  1.57e-4  then 
mcolor  -  2 

else 

mcolor  -  1 
end  if 

pset  (x/mitokm,y/fttom) ,  mcolor 

else 

if  mgrad  <  0  then  pset  (x/mitokm, y/f ttom) ,  4 
end  if 

x  -  x  +  xinc 

if  x  >  of f set ( i+1 )  then  i  -  i  +  1 

wend 

y  —  y  +  yinc 
i  -  1 

if  y  >  p(i,n+l,l)  then  n  -  n  +  1 
j  “  n 

wend 

locate  25,1 

if  mode$  O  "EGA16"  then  print  " 


'  put  the  dots  In  an  array  for  later  display 

window 

view 

get  (12. 5*charwd+0. l*charwd,  2. 4*charht+0. l*charht)_ 

-  (xcol*charwd-0. l*charwd,  yrow*charht-0. l*charht) ,  dots 

else 


window 

view 

put  (12 . 5*charwd+  0.  l*charwd,  2. 4*charht+0. l*charht) ,  dots,  xor 
end  if 

'  establish  viewport  and  window  for  drawing 

view  (12.5  *  charwd,  2.4  *  charht)  -  (xcol  *  charwd,  yrow  *  charht) 
window  (O.hmax)  -  (xmax.O) 

'  find  the  layer  containing  the  transmitter 

ok$  -  "no" 
i  -  1 

while  ok$-  "no  "  and  i  <-  nmax(l) 
if  tran  <  p(l,i,l)  then 
i  -  i  -  1 
ok$  -  "yes" 

elseif  tran  -  p(l,i,l)  then 
ok$  -  "yes" 

else 

i  -  i  +  1 
end  if 

wend 

tlevel  -  i 

for  j  -  1  to  nrays 

'  initialize  the  beginning  values 

n  -  1 

i  -  tlevel 

bht  -  tran 

eht  -  bht 

brng  -  0 

erng  -  0 

rngnow  -  rlim(n) 

rnginc  -  rlim(n)  /  25 

if  j  -  1  then 

balpha  -  startangle 

else 

balpha  -  startangle  +  (j  -1)  *  angleinc 
end  if 


loop  until  maximum  height  or  range  is  reached 


while  eht  <  hmax  /  fttom  and  erng  <  xlimit  and  i  <  nmax(l) 

compute  local  M-unit  gradient,  ending  range  and  angle 

mgrad  -  cm(n,i)  *  (brng-offset (n) )  +  p(n,i,3) 
if  (mgrad  -  0)  then  mgrad  -  l.e-6 
erng  -  brng  +  rnginc 

if  (erng  >  rngnow)  then  erng  -  rngnow 
ealpha  -  balpha  +  mgrad  *  (erng  -  brng) 

check  to  see  if  ray  has  passed  through  a  maximum  or  minimum 
if  so,  compute  a  new  ending  range  and  angle 

if  (balpha<0  and  ealpha>-0)  or  (balpha>0  and  ealpha<=0)  then 
ealpha  -  0 

erng  -  brng  -  balpha  /  mgrad 
end  if 

compute  an  ending  height  and  the  boundary  heights  at  this 
ending  range 


eht  -  bht  +  (ealphaA2  -  balphaA2)  /  (2.e-3  *  mgrad) 
topht  -  ch(n, i+1)  *  (erng-offset(n) )  +  p(n,i+l,l) 
btmht  -  ch(n, i)  *  (erng-of fset(n) )  +  p(n,i,l) 

if  ray  has  penetrated  layer  boundary,  compute  range,  height,  and 
angle  of  boundary  crossing 

if  eht  >  topht  or  eht  <  btmht  then 
if  eht  >  topht  then 
11  -  i  +  1 
flag  -  1 

else 

11  -  i 
flag  -  0 
end  if 

a  -  mgrad  /  2  .  e  -  3 
b  -  2  *  balpha  /  2.e-3  -  ch(n,ll) 

c  -  bht  -  ch(n,ll)  *  (brng  -  offset(n))  -  p(n,ll,l) 
delr  -  (-b  +  sqr(bA2  -  U  *  a  *  c))  /  (2  *  a) 

if  delr  <  0  or  delr  >  rnginc  then_ 

delr  -  (-b  -  sqr(b  2  -  U  *  a  *  c))  /  (2  *  a) 
if  delr  -  0  then  delr  -  rnginc 
erng  -  brng  +  delr 

ealpha  -  balpha  +  mgrad  *  (erng  -  brng) 
eht  -  ch(n,ll)  *  (erng  -  offset(n))  +  p(n,ll,l) 
if  flag  -  1  then 
i  -  i  +  1 

else 

i  -  i  -  1 
end  if 
end  if 
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reverse  angle  if  ray  has  reached  the  ground 


if  (eht  <-  0)  then 
i  -  1 

ealpha  -  -ealpha 
eht  -  0 
end  if 

'  draw  the  ray  segment  on  the  screen 

line  (brng/mitokm,  bht/fttom)  -  (erng/mitokm,  eht/f ttom) , 14 

'  check  to  see  if  ray  has  reached  next  profile  range.  If  so, 

'  increment  range  counter  and  get  new  step  size 

if  erng  -  rngnow  and  erng  <  xlimit  then 
n  -  n  +  1 

rngnow  -  rngnow  +  rlim(n) 
rnginc  -  rlim(n)  /  25 
end  if 

'  reinitialize  starting  height,  angle,  and  range 

bht  -  eht 
balpha  -  ealpha 
brng  -  erng 

wend 
next  j 

'  toggle  background 


window 

view 

locate  25,1 

print  "Press  <t>  to  toggle  background  or  <ENTER>  to  continue."; 

toggle: 
a$  -  "" 
while  a$  -  "" 
a$  -  inkey$ 

wend 

if  a$  -  chr$(84)  or  a$  -  chr$(116)  then_ 

put  (12 . 5*charwd+  0.  l*charwd,  2. 4*charht+0.  l*charht)  ,dots,  xor 
if  a$  -  chr$(13)  then  goto  continue 
goto  toggle 


continue : 


Subroutine:  gradient 


'  Purpose:  1.  Compute  an  M-unit  gradient 

sub  gradient  static 
for  i  -  1  to  nprof 

'  compute  the  M-unit  value  at  one  kilometer 

ok$  -  "no" 
j  -  nmax( i) 

if  p(i , j ,1)  <  1000  then 

diff  -  (1000  -  p(i , j , 1) )  /  1000 
matkm  -  (118  *  diff)  +  p(i,j,2) 
elseif  p(i,j,l)  -  1000  then 
matkm  -  p( i , j , 2) 

else 

while  ok$  -  "no" 

if  p(i,j,l)  >  1000  then 
J  -  j  -  1 

else 

matkm  -  p(i,j+l,2)  -  (p(i,j+l,2)  -  p(i  ,  j , 2) )  *_ 

(pd.J+l.l)  -  1000)  /  (p(i,j+l,l)  -  p(  i  .  j  .  1) ) 
ok$  -  "yes" 
end  if 

wend 
end  if 

'  compute  the  M-unit  gradient  for  each  profile 

for  j  -  1  to  nmax(i)  -  1 

delu  -  p(i , j+1 , 2)  -  p(i , j ,2) 
delh  -  p(i , j+1 , 1)  -  p(i,j,l) 
if  delh  O  0  then 

p(i , j , 3)  -  l.e-3  *  delu  /  delh 

else 

P(i  .  j  .3)  -  l.e-6 
end  if 
next  j 

'  compute  the  M-unit  gradient  at  the  top  of  the  profile 

matsfc  -  p ( i , 1,2) 

h  -  p(i ,nmax( i) , 1)  /  1000 

a  -  log(matsfc  /  (matkm  -  157)) 

p( i , nmax( i) , 3)  -  157  -  (matsfc  *  a  *  exp(-a  *  h)) 
next  i 
end  sub 
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